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ABSTRACT 

Accurate  modeling  of  Space  Surveillance  sensors  is  necessary  for  a  variety  of  applications.  Accurate  models  can  be 
used  to  perform  trade  studies  on  sensor  designs,  locations,  and  scheduling.  In  addition,  they  can  be  used  to  predict 
system-level  performance  of  the  Space  Surveillance  Network  to  a  collision  or  satellite  break-up  event. 

A  high  fidelity  physics-based  radar  simulator  has  been  developed  for  Space  Surveillance  applications.  This 
simulator  is  designed  in  a  modular  fashion,  where  each  module  describes  a  particular  physical  process  or  radar 
function  (radio  wave  propagation  &  scattering,  waveform  generation,  noise  sources,  etc.)  involved  in  simulating  the 
radar  and  its  environment.  For  each  of  these  modules,  multiple  versions  are  available  in  order  to  meet  the  end-users 
needs  and  requirements.  For  instance,  the  radar  simulator  supports  different  atmospheric  models  in  order  to  facilitate 
different  methods  of  simulating  refraction  of  the  radar  beam.  The  radar  model  also  has  the  capability  to  use  highly 
accurate  radar  cross  sections  generated  by  the  method  of  moments,  accelerated  by  the  fast  multipole  method.  To 
accelerate  this  computationally  expensive  model,  it  is  parallelized  using  MPI. 

As  a  testing  framework  for  the  radar  model,  it  is  incorporated  into  the  Testbed  Environment  for  Space  Situational 
Awareness  (TESSA).  TESSA  is  based  on  a  flexible,  scalable  architecture,  designed  to  exploit  high-performance 
computing  resources  and  allow  physics-based  simulation  of  the  SSA  enterprise.  In  addition  to  the  radar  models, 
TESSA  includes  hydrodynamic  models  of  satellite  intercept  and  debris  generation,  orbital  propagation  algorithms, 
optical  brightness  calculations,  optical  system  models,  object  detection  algorithms,  orbit  determination  algorithms, 
simulation  analysis  and  visualization  tools.  Within  this  framework,  observations  and  tracks  generated  by  the  new 
radar  model  are  compared  to  results  from  a  phenomenological  radar  model.  In  particular,  the  new  model  will  be 
used  to  simulate  an  S-band  upgrade  to  the  space  fence. 

1.  INTRODUCTION 

Computer  models  of  radars  in  the  Space  Surveillance  Network  can  be  a  useful  tool  for  evaluating  new  sensor 
locations,  designs,  and  operational  parameters.  Such  models  can  also  be  used  to  study  tasking  algorithms  and 
predict  the  system  performance  to  new  debris.  However,  when  performing  many  trade  studies  or  simulations 
spanning  several  simulated  days,  these  models  can  have  unacceptable  run  times.  This  leads  to  a  trade-off  between 
the  inclusion  of  all  relevant  physics  in  the  model  and  simulation  time.  This  paper  details  a  radar  model  which 
answers  the  run-time  problem  while  maintaining  a  high  level  of  modeling  fidelity  by  using  parallelization.  The 
radar  model,  which  is  implemented  in  Lawrence  Livermore  National  Laboratory’s  TESSA  (Testbed  Environment 
for  Space  Situational  Awareness)  is  designed  to  run  on  Linux  cluster  supercomputers.  It  utilizes  a  modular 
framework  which  allows  different  physics  models  to  be  used  as  appropriate  for  any  radar.  This  general  nature 
makes  it  a  good  tool  for  modeling  many  different  sensors  with  the  same  software.  Highly  accurate  modules  have 
been  created  for  noise  sources,  antenna  effects,  atmospheric  effects  on  propagation,  and  the  radar  cross  sections  of 
targets.  These  modules  are  combined  in  a  framework  which  uses  parallelization  to  divide  the  simulation  work  load 
on  a  pulse-by-pulse  basis  with  each  processor  being  responsible  for  all  computations  involving  a  single  pulse.  This 
paper  will  first  describe  several  of  the  physics  modules  in  the  simulation,  then  explain  the  parallelization  scheme  and 
present  an  example  simulation  with  timing  results. 

2.  RADAR  MODELING 
2.1  EXTERNAL  NOISE  SOURCES 

One  of  the  fundamental  quantities  that  a  radar  can  measure  with  its  antenna  is  the  radiation  power  per  unit  area,  per 
unit  frequency.  In  the  microwave  regime,  this  is  usually  given  the  name  brightness  temperature,  Tp .  There  are 
various  sources  which  contribute  to  the  brightness  temperature.  They  include  external  sources  outside  of  the  earth’s 
atmosphere  (e.g.  cosmic  or  galactic  noise),  emission  from  atmospheric  gases,  and  emission  and  scattering  from  the 
ground.  This  section  will  deal  with  modeling  cosmic  contributions  to  the  brightness  temperature, 
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—  Tcms  -h  7^[ .  Tcmb  is  the  temperature  due  to  the  cosmic  microwave  background  emission  and  is  usually 
modeled  as  a  constant,  2.73  K.  A  common  way  to  model  the  galactic  noise,  7^,  is  to  use  a  well  defined  median 
temperature  at  a  reference  frequency,  and  apply  a  scaling  law  to  model  the  galactic  noise  temperature  at  a  different 
frequency,  as  is  shown  in  eqn.  2.1.1. 
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Typical  values  for  the  parameters  in  equation  2.1.1  are:  Tr^r-  =20  K,  ft  =  2.75  ,  and  fr€f  =  408  MHz.  This  method 
is  very  computationally  efficient,  but  effectively  averages  over  any  variations  in  temperature  due  to  galactic 
structure,  bright  point  sources,  etc. . .  In  [1],  a  method  for  computing  the  sky  temperature  is  provided  that  models  the 
angular  distribution  of  the  galactic  noise  sources  and  covers  a  wide  range  of  frequencies  from  10  MHz  to  300  GHz. 
The  model  is  based  on  a  numerical  fit  to  a  radio  emissions  data  set  derived  from  all  publically  available  total  power, 
large-area  radio  surveys.  The  inputs  to  the  galactic  nosie  model  include  the  radar’s  operating  frequency,  location, 
year,  month  and  day,  and  time.  The  outputs  are  contained  in  a  table  of  temperatures  indexed  with  galactic  latitude 
and  longitde.  A  routine  to  convert  from  galactic  longitude  and  latitude  is  used  to  translate  the  results  to  the  radar’s 
coordinate  system.  Using  this  method  allows  a  detailed  model  of  the  angular  distribution  of  the  galactic  noise  as 
measured  by  a  particular  radar  at  any  location,  date,  time  of  day  and  operating  frequency.  An  example  of  this  is 
shown  below  in  Figure  1,  for  a  radar  located  at  Death  Valley,  CA  at  night  time,  operating  at  a  frequency  of  408 
MHz.  Galactic  radio  emissions  can  be  a  significant  source  of  noise  for  frequencies  below  ~1  GHz,  but  rapidily 
become  neglibible  for  frequencies  much  beyond  this  (see  fig.  1,  right  panel).  Future  work  in  this  area  will  focus  on 
modeling  man-made  sources  of  RF  noise  and  scattering  from  the  ground. 
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Figure  1  (Left  Panel)  Simulated  galactic  noise  distribution  at  408  MHz  for  a  location  in  Death  Valley.  (Right 
Panel)  The  galactic  noise  temp,  dependence  with  frequency,  taken  over  all  gal.  latitudes  and  longitudes. 


2.2  ATMOSPHERIC  MODELS 


This  section  will  describe  the  different  atmosphere  models  that  are  currently  available  in  the  simulation  framework. 
The  Earth’s  atmosphere  can  affect  the  propagation  of  Radio  Frequency  (RF)  energy  in  a  number  of  different  ways, 
e.g.  refractive  bending,  attenuation,  range-delay,  Doppler  frequency  shifts,  polarization  rotation,  frequency 
spreading  and  rapid  amplitude  and  phase  scintillation  [2] .  Accurately  modeling  these  effects  is  a  critical  step  in  the 
simulation,  as  it  affects  the  quality  of  the  measurements  generated  by  the  radar  which  are  used  in  computing  orbital 
parameters  of  objects  being  tracked  by  the  radar.  In  order  to  support  simulation  activities  that  require  a  high  level  of 
accuracy  in  the  modeling,  we  have  incorporated  high  fidelity  models  of  the  atmosphere  into  the  simulation  (when 
available).  Lower  fidelity  models  are  also  included  in  the  simulation,  which  are  useful  for  conducting  a  large  number 
of  simulations  very  quickly.  The  inputs  to  the  different  models  can  include  geographical  latitude,  longitude,  altitude, 
time  and  date.  The  outputs  for  the  neutral  atmosphere  models  are  altitude  dependent  temperature,  pressure  and  water 
vapor  pressure.  The  output  for  the  ionosphere  model  is  the  electron  density  as  a  function  of  the  model  parameters. 
The  simulated  atmosphere  is  then  used  in  modeling  the  effect  of  the  atmosphere  on  the  propagation  of  the  RF 
energy. 

Currently,  there  are  7  neutral  atmospheric  models  and  2  ionospheric  models  available  for  use  within  the  simulation 
framework.  The  first  type  of  neutral  atmosphere  model  available  to  the  simulation  user  is  an  exponential 
approximation  to  the  neutral  atmosphere,  see  [3]  for  details.  The  second  class  of  neutral  atmosphere  model  includes 
a  set  of  Standard  Atmospheres  [4].  A  “Standard  Atmosphere”  is  defined  as  a  vertical  distribution  of  atmospheric 


parameters  that  are  representative  of  Earth’s  atmosphere.  There  are  6  types  of  Standard  Atmospheres  in  the 
simulation  framework;  US  Standard,  Tropical,  Mid-latitude  (summer  and  winter),  and  Sub- Arctic  (summer  and 
winter).  Each  standard  atmosphere  makes  use  of  a  data  base  that  provides  realistic  vertical  profiles  of  atmospheric 
temperature,  pressure,  water  content  and  gas  mixing  ratios  from  ground  level  to  an  altitude  of  120  km.  These 
parameters  are  used  in  computing  the  index  of  refraction  for  modeling  propagation  effects  in  the  neutral  atmosphere; 
see  section  2.3  for  more  details. 

Two  different  models  of  the  ionosphere  are  available  to  the  simulation  user.  The  first  [5]  approach  is  a 
phenomenological  model  of  the  ionosphere  that  captures  large  scale  variations  in  the  electron  density  as  a  function 
of  altitude.  The  second  ionosphere  model  available  to  the  user  is  the  International  Reference  Ionosphere  (IRI)  [6]. 
The  IRI  modeling  effort  was  started  in  1969  and  has  been  steadily  improved  on  each  year.  This  provides  a  higher 
fidelity  estimate  of  the  electron  density  as  a  function  of  a  large  number  of  variables.  Typical  runs  will  execute  the 
IRI  routines  using  the  radar’s  position,  simulation  date  and  the  time  of  day  with  all  other  IRI  variables  to  their 
default  values.  An  example  of  the  IRI  model  output  is  shown  in  Figure  2.  In  the  left  panel  the  Total  Electron  Content 
(TEC)  is  plotted  as  a  function  of  latitude  and  longitude.  In  the  right  panel,  range  delays  due  to  the  ionosphere  are 
plotted  as  a  function  of  frequency.  For  operating  frequencies  above  a  few  GHz,  ionospheric  effects  are  typically  not 
important  and  can  be  readily  corrected  for  by  the  radar.  For  operating  frequencies  much  below  1  GHz,  ionospheric 
effects  become  more  pronounced  and  will  require  the  radar  to  perform  a  dual  frequency  correction  or  employ 
sophisticated  single  frequency  correction  schemes.  These  sorts  of  corrections  performed  by  the  radar  are  not 
currently  modeled  in  the  simulation  and  will  be  a  focus  for  future  development  work.  In  addition  to  modeling  the 
corrections  performed  by  the  radar,  future  work  will  focus  on  incorporating  Numerical  Weather  Models  into  the 
simulation  framework.  This  will  enable  a  more  realistic  model  of  the  day-to-day  (and  even  hour-to-hour)  changes  in 
the  atmospheric  conditions  that  affect  the  operations  of  the  radar. 


2002  April  25  tecu  Range  Delay  Trends 


Figure  2  (Right  panel)  Ionospheric  TEC  as  a  function  of  latitude  and  longitude  for  April  25,  2002,  UT  =  5 
hours.  (Left  panel)  Ionospheric  range  delay  as  a  function  of  frequency  for  a  number  of  different  conditions. 


2.3  RF  PROPAGATION  EFFECTS 

The  atmospheric  conditions  are  simulated  for  a  given  radar’s  location,  current  time,  day,  month  and  year  using  one 
of  the  models  described  in  section  2.2.  The  temperature,  pressure,  water  vapor  pressure  and  electron  density 
obtained  from  the  atmosphere  model  are  used  to  compute  a  3  dimensional,  complex  refractivity  map  of  the  sky  as 
viewed  from  the  radar’s  position.  This  refractivity  map  is  then  used  with  a  ray-tracing  procedure  to  compute  the  path 
taken  by  the  RF  ray  as  it  travels  through  the  atmosphere.  Once  the  path  has  been  computed,  the  amount  of  refraction 
that  occurred  is  known,  the  range  delay  can  be  computed  and  the  total  attenuation  suffered  by  the  ray  can  be 
calculated  by  integrating  the  specific  attenuation  along  the  ray’s  path.  The  atmospheric  attenuation  will  affect  the 
measured  Radar  Cross  Section  (RCS)  and  the  signal-to-noise  ratio.  The  range  delay  and  refractive  bending  will 
affect  the  measurement’s  range  and  angular  accuracies. 

Liebe’s  Millimeter- wave  Propagation  Model  (MPM)  [7]  is  used  to  model  the  complex  refractivity  in  the  neutral 
atmosphere  for  frequencies  up  to  1000  GHz.  The  complex  refractivity,  A^C/),  as  defined  by  Leibe  is  composed  of 
three  parts,  a  nondispersive  component  which  is  real  and  positive,  along  with  a  set  of  real  Nip  'and  imaginary 
N  if)  ' '  components, 

NXf}t  =NQ+N(pr-jN(p”  ppm. 

The  imaginary  component  is  related  to  power  attenuation  through  the  following  expression, 

a  =  0.1S2G  fNTT(p  dB/km 


and  the  real  part  describe  phase  dispersion, 

fj  =  1,2008  deg/km , 

where  f  is  the  frequency  in  GHz.  An  example  of  how  the  index  of  refraction  is  used  in  the  code  base  is  shown  in 
figure  3.  The  left  panel  show  a  comparison  of  the  refractive  bending  in  the  neutral  atmosphere  as  a  function  of  the 
initial  elevation  angle.  Significant  bending  occurs  for  low  elevations  and  moist  atmospheres  (Tropical).  The  middle 
panel  shows  the  index  of  refraction  in  the  neutral  atmosphere  as  a  function  of  altitude.  Variations  in  the  refractive 
index  due  to  different  atmospheric  conditions  are  limited  to  the  first  ~30  kilometers  of  the  earth’s  surface.  The  ray 
tracing  routine  was  validated  using  a  number  of  different  cases.  The  left  panel  shows  the  refractive  bending  in  an 
exponential  atmosphere  as  computed  using  the  ray  trace  routine  (black  asterisk)  and  an  analytical  approximation 
(black  line)  [8],  showing  excellent  agreement.  Another  example  of  the  ray  trace  validation  is  show  in  the  rightmost 
panel.  A  simulation  of  a  Maxwell’s  Eye  refractive  medium  was  performed  (grayscale  indicates  index  of  refraction) 
and  a  number  of  different  rays  were  traced.  The  comparison  between  the  numerical  results  and  theoretical 
expectations  agreed  to  within  0.001%  for  this  simple  case. 
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Figure  3.  (Top  Left)  A  comparison  of  the  elevation  dependent  refractive  bending  for  a  number  of  different 
cases.  (Top  Right)  Index  of  refraction  as  a  function  of  altitude.  (Bottom)  An  example  of  a  Ray  Trace 
validation  case. 


The  Appleton-Hartree  equation  (see  [9])  serves  as  the  starting  point  for  modeling  the  index  of  refraction  in  the 
ionosphere, 

_  V 

n*=  1  - 


Y2s\n29  ,  {Y*s in*e  ■,  nV'2 


where  a  =  oj^f  &j2  ,  Y  =  sfmw,  6T  =  1  —  and  s2/  £^Sn.  The  angular  plasma  frequency  is  , 


N  is  the  electron  density,  e  and  m  are  the  charge  and  mass  of  an  electron,  jUd  are  the  electric  and  magnetic 
permittivity  of  free  space,  H  is  the  geomagnetic  field  intensity,  v  is  the  collision  frequency,  0  is  the  angle  made 
between  the  wave  normal  and  the  geomagnetic  field,  and  the  +  refers  to  the  ordinary  and  extra-ordinary  waves.  As 
described  in  [9],  a  simplification  can  be  made  for  frequencies  above  VHF,  yielding  the  expression 

n  ss  1--Y  =  1 

2  tLl‘ 

where  b  =  /2  E$m.  This  approximation  ignores  the  effect  of  the  splitting  of  the  ordinary  and  extra-ordinary  waves 

and  thus  reduces  the  index  of  refraction  to  a  function  solely  based  on  the  electron  density,  radar  operating  frequency, 
and  a  few  physical  constants.  The  electron  density  can  be  simulated  using  the  ionosphere  models  described  in  sec. 
2.2.2  and  will  be  a  function  of  year,  month,  day,  time  of  day  ,  geographic  latitude,  longitude  and  altitude 


Since  the  earth’s  atmosphere  emits  and  absorbs  RF  energy,  we  use  Radiative  Transfer  Theory  to  model  the  noise 
power  measured  with  the  radar  antenna,  Tb.  Following  [10],  the  noise  power  can  be  related  to  the  galactic  brightness 
temperature  Tgal,  and  the  atmospheric  temperature  T. 

rb  =  Tgci  a  + ;0L  «(xr)nxr)3-<°^dx', 

where  ir(xr.,  x)  =  \x,  h(sD  ds  is  the  optical  thickness  and  T(xT ]  is  the  temperature  of  the  atmosphere  at  an  altitude  x' . 
The  ray  trace  routine  is  used  within  the  integration  routine.  By  ignoring  scattering  effects  we  limit  ourselves  to  cases 
without  significant  rain,  smoke,  and  other  scattering  sources.  An  example  of  the  effects  of  radiative  transfer  theory 
has  on  the  noise  power  measured  by  the  radar  is  shown  in  figure  4  below. 


Frequency  (GHz) 

Figure  4.  Example  of  a  Radiative  Transfer  Theory  calculation  for  a  number  of  different  elevation  angles. 


2.4  ANTENNA  EFFECTS 

The  radar  antenna  is  used  to  transmit  and  receive  RF  energy,  and  can  have  many  different  designs  and  capabilities. 
Currently,  the  simulation  supports  rectangular,  phased  array  antennas.  The  option  of  using  a  dish  antenna  will  be 
incorporated  into  the  simulation  as  part  of  future  work.  The  phased  array  antenna  is  defined  in  the  simulation  by  a 
number  of  variables;  (1)  type  of  radiator,  (2)the  number  of  radiators  along  the  width  and  length  of  the  antenna,  (3) 
the  spacing  of  the  radiators,  (4)  orientation  of  radiators,  (5)  the  number  of  bits  used  in  the  phase  shifters,  (6)  the 
location  of  the  antenna,  (7)  orientation  of  antenna,  (8)  antenna  efficiency,  and  (9)  the  drift  rate  of  Inertial 
Measurement  Unit  (IMU). 

The  user  can  select  from  two  pre-defined  radiating  elements  (half-wave,  Hertzian  dipole).  If  desired,  the  option  of  a 
user-defined  radiating  element  (such  as  a  waveguide,  horn,  etc. . .)  is  supported  by  the  simulation  as  well.  The  user 
has  the  option  of  switching  on  a  model  of  a  strap-down  IMU  co-located  with  the  radar  antenna.  This  will  allow  a 


simulation  of  the  performance  degradations  due  to  imperfect  motion  compensation,  errors  in  registration,  etc. . .  The 
co-  and  cross-polarization  gains  of  the  antenna  are  computed  as  a  function  of  the  antenna’s  boresite  and  scan  angle. 
Following  the  discussion  in  [10],  the  noise  power  measured  by  the  antenna  is  computed  using  the  following  formula, 

PJJ.f;9^\rQ)  sin  9  d9dd> 

9  d9dd 

where  Tp  is  defined  in  section  2.3  and  Pn  is  the  total  antenna  pattern  (co-  and  cross-polar).  The  calculation  of  ^nfjis 
computationally  intensive  so  it  is  split  into  two  parts.  Tp  varies  slowly  with  time  and  can  be  expensive  to  compute, 
requiring  several  minutes  if  a  raytracing  radiative  transfer  method  is  used.  Because  it  is  slow  to  compute,  a  lookup 
table  of  Fp  is  created  over  the  whole  visible  sky  at  one  degree  intervals.  This  table  filled  once  and  used  for  up  to  one 
hour  of  simulation  time  before  being  recomputed.  The  antenna  gain  Fn  (/,  0 ,  0  jfD)  changes  at  every  pulse  as  the 
antenna  scans  the  sky,  but  it  can  be  computed  quickly  for  each  pulse,  then  integrated  against  the  precomputed  table 
of  sky  temperatures. 


tSf  F:  )  — 


2.5  RADAR  CROSS  SECTION  CALCULATION 

Having  an  accurate  radar  cross  section  model  for  targets  is  crucial  to  performing  an  accurate  radar  simulation. 
TESSA  has  several  radar  cross  section  models  which  can  be  used  depending  on  the  desired  accuracy  and  availability 
of  information.  When  catalog  objects  are  used  in  the  simulation,  their  RCS  value  from  the  Satellite  Situation 
Report  is  used  if  there  is  no  mesh  information  available.  For  other  objects,  either  satellites  with  CAD  models 
available  or  simulated  debris  objects  from  a  hydro  code,  an  RCS  calculation  method  is  chosen  based  on  the  objects’ 
size.  For  very  small  debris  objects  where  the  object  size  is  much  less  than  a  wavelength,  a  Rayleigh  scattering 
approximation  is  used  [11].  For  larger  objects  in  the  resonant  regime,  EIGER,  Electromagnetic  Interactions 
GEneRalized,  an  in-house,  parallel,  method-of-moments  solver  is  used  [12].  This  solver  supports  metals  and 
dielectric  materials,  and  provides  highly  accurate  solutions  for  RCS.  For  objects  which  are  many  wavelengths  in 
size,  two  models  can  be  used.  One  is  a  code  based  on  Physical  Optics,  which  reduces  the  current  integrals  over 
facets  to  integrals  over  their  edges  [13]  along  with  the  Physical  Theory  of  Diffraction  [11]  to  handle  scattering  from 
edges.  The  other  computation  method  used  for  large  scatterers  is  the  fast  method  implementation  of  EIGER. 

EIGER  contains  a  parallel  implementation  of  the  Fast  Multipole  Method  -  Fast  Fourier  Transform  (FMM-FFT) 
method  [14].  This  method  allows  for  highly  accurate  fast  solutions  for  radar  cross  section  problems,  while  making 
load  balancing  and  communication  in  parallel  simpler  than  the  Multilevel  Fast  Multipole  Method  [15].  Figure  5 
shows  a  comparison  of  the  RCS  generated  using  the  Physical  Optics  method  and  EIGER.  While  the  Physical  Optics 
has  good  agreement  with  the  high  fidelity  model  for  specular  scattering,  returns  near  monostatic  vary  by  lOdB 
between  the  methods.  It  is  important  to  note  that  the  difference  in  computational  cost  between  the  methods;  EIGER 
required  1.7  hours  on  16  processors  for  the  simulation,  while  the  Physical  Optics  method  took  just  under  10  minutes 
on  a  single  processor. 

For  RCS  computations  using  EIGER  or  Physical  Optics  a  cumulative  distribution  of  RCS  values  for  each  object  is 
precomputed  before  the  radar  simulation.  This  distribution  is  formed  by  computing  the  monostatic  RCS  of  each 
object  at  1536  equally-spaced  locations  and  filling  the  results  into  a  histogram.  During  the  simulation,  when  the 
RCS  of  a  target  is  needed,  it  is  taken  randomly  from  the  cumulative  distribution.  Radar  Cross  Section  values  which 
are  generated  using  Rayleigh  scattering  or  from  the  Satellite  Situation  Report  are  not  formed  into  a  distribution; 
instead  a  single  value  is  always  used  for  the  object. 


BistaticRCS  for  Cosmos  2251  at  5  GHz 
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Figure  5.  The  Bistatic  RCS  computed  in  EIGER  and  with  Physical  Optics  for  the  Cosmos  2251  satellite  at  5 
GHz  and  an  incident  angle  of  150  degrees.  The  two  methods  show  reasonable  agreement  for  the  main 
specular  peeks,  but  not  for  the  backscattered  result  near  150  degrees.  The  incident  field  direction  is  shown  by 
the  red  arrow. 


2.6  MPI  PARALLELIZATION 

The  computationally  intensive  nature  of  this  radar  model  demands  parallelization.  By  parallelizing  the  model,  we 
can  decrease  the  run  time  so  that  trade  studies  of  radar  designs  or  placement  can  be  run  quickly  without  sacrificing 
simulation  accuracy.  The  radar  simulation  is  parallelized  in  several  stages.  All  parallelization  is  done  using  the 
Message  Passing  Interface  (MPI)  [16].  The  TESSA  framework  that  this  model  is  a  part  of  automatically  assigns 
each  radar  in  a  simulation  to  its  own  processor  or  bank  of  processors.  If  low  fidelity,  inexpensive  modules  alone  are 
used,  this  level  of  parallelization  may  be  sufficient.  When  more  intensive  modules  are  used,  each  radar  is 
parallelized  further.  First,  the  simulation  duration  is  broken  into  smaller  time  blocks,  typically  one  hour  in  duration. 
For  these  short  blocks,  some  of  the  modules  fill  lookup  tables  which  are  referenced  at  each  pulse  rather  than  being 
recomputed.  Examples  of  this  include  the  ionosphere  module,  which  samples  the  electron  density  with  5km 
spacing,  and  the  antenna  module,  which  performs  part  of  the  antenna  temperature  calculation  through  a  table. 

These  tables  can  require  many  minutes  to  fill,  and  the  same  data  is  needed  on  every  processor.  To  handle  these 
tables  efficiently,  these  tables  are  split  with  each  processor  computing  the  same  number  of  entries.  The  completed 
entries  are  then  shared  so  a  complete  table  is  on  every  processor. 

After  the  tables  are  completed,  the  radar  pulses  are  split  among  the  processors.  Each  processor  is  responsible  for  the 
complete  simulation  of  a  contiguous  block  of  pulses.  Dividing  the  workload  by  pulses  has  two  main  advantages. 

The  first  is  that  most  of  the  individual  modules  that  are  called  at  each  pulse  do  not  have  to  be  parallelized,  which 
speeds  development.  The  second  advantage  is  for  saving  time  on  orbital  propagation.  Orbital  propagation  of  radar 
targets  can  be  time  consuming  if  a  force  model  propagator  is  used.  The  force  model  propagator  in  TESSA 
precomputes  a  table  of  locations  at  five  to  ten  minute  intervals,  then  propagates  the  objects  only  from  the  last  table 
entry  during  the  simulation.  By  splitting  the  pulses  in  time,  each  processor  only  has  to  update  the  position  and 
process  from  the  nearest  table  entry. 

Once  all  pulses  are  simulated,  the  observations  are  packed  and  sent  to  the  master  processor  for  assignment  to  tracks. 
While  this  final  stage  is  not  load  balanced,  the  simple  tracking  algorithm  used  is  not  a  dominant  cost  in  the 
simulation.  Using  this  method,  efficient  scaling  has  been  demonstrated  from  1-16  processors  per  radar.  Future 
work  in  parallelization  will  include  parallelizing  advanced  tracking  algorithms  as  they  are  implemented. 

3.  SPACE  FENCE  EXAMPLE 

To  demonstrate  the  models  capabilities,  it  was  used  to  implement  a  notional  S-band  space  fence  radar  located  at  four 
installations.  This  radar  is  not  meant  to  represent  a  specific  future  space  fence  radar,  but  to  be  similar  to  one 
possible  design.  The  radar  sites  included  three  locations  in  the  United  States  roughly  along  the  33rd  parallel  (Gila 
River,  Lake  Kickapoo,  and  Jordan  Lake)  as  well  as  one  Australian  site,  Canbera,  located  at  latitude  35.3  South.  The 
radar  operating  frequency  was  3.5  GHz,  with  an  opening  fan  in  the  East-West  direction  from  10  degrees  above 


either  horizon  and  the  ability  to  track  objects  through  22.5  degrees  of  arc  when  passing  through  the  fan.  The  radar  is 
assumed  to  use  a  monopulse  array  for  the  actual  tracking  of  objects  once  detected.  The  radar  is  designed  to  detect 
objects  as  small  as  5cm  at  a  range  of  1000  km,  and  has  a  maximum  effective  range  of  20,000  nautical  miles.  These 
parameters  give  it  better  sensitivity  and  slightly  better  range  than  the  existing  VHF  space  fence,  while  maintaining  a 
similar  level  of  accuracy.  Because  this  radar  model  operates  in  monostatic  mode  with  some  ability  to  track  an  object 
that  penetrates  its  beam,  preliminary  object  orbits  can  be  determined  on  the  first  orbit. 

Detection  and  accuracy  of  a  given  observation  depends  on  the  signal -to-noise  ratio  of  the  observation.  The  radar  is 
tuned  to  have  a  skin  track  range  of  5225  km,  which  allows  it  to  detect  a  20cm  object  at  2200km  or  a  5cm  object  at 
around  1000km.  The  minimum  signal-to-noise  ratio  used  for  detection  is  12  dB.  During  the  simulation,  the  signal  - 
to-noise  ratio  for  a  given  target  is  used  to  determine  the  error  in  the  observation,  which  is  assumed  to  be  Gaussian. 

A  Cramer-Rao  lower  bound  model  for  error  from  monopulse  radar  is  used.  This  model  results  in  angular  accuracy 
of  around  l/6th  of  the  beam  width  for  12  dB  targets  at  the  limit  of  detection,  and  around  l/30th  of  a  beamwidth  for 
extremely  strong  targets. 

During  the  simulation,  the  radar  is  scheduled  to  continuously  observe  the  sky.  A  simple  pulse-scheduling 
algorithm  was  assumed,  which  allowed  the  radar  to  update  the  position  of  every  object  within  its  field  of  view  once 
a  second.  For  every  observation,  the  signal-to-noise  ratio  was  computed  taking  into  account  radar  parameters,  such 
as  power  and  antenna  pattern,  distance  to  target,  and  object  radar  cross  section.  This  signal-to-noise  ratio  then 
determined  the  accuracy  of  the  observation  and  if  the  object  was  detected.  For  all  objects  detected,  the  observation 
plus  appropriate  error  was  passed  to  the  tracking  algorithm. 

The  test  runs  were  performed  on  HERA,  a  LINUX -based  supercomputer  cluster  at  Lawrence  Livermore  National 
Laboratory.  HERA  consists  of  13,824  2.3  GHz  AMD  Opteron  Quad  Core  processors  with  2  GB  of  memory  per 
processor.  The  simulation  required  only  a  small  fraction  of  the  computer,  using  between  4  and  64  processors. 

The  simulated  space  fence  radars  were  examined  for  their  ability  to  detect  debris  from  the  Cosmos  225 1  -  Iridium 
33  collision  on  February  10,  2009.  This  debris  was  created  using  a  hydrodynamic  code  and  full  CAD  models  of  the 
two  satellites.  The  debris  catalog  contained  414  objects  with  mean  RCS  values  ranging  from  3.7  dbsm  to  -57  dbsm. 
The  radar  simulation  was  carried  out  at  two  levels  of  fidelity.  Both  used  the  antenna  model  and  high  accuracy  RCS 
values.  However,  for  case  one,  all  galactic  noise  sources,  atmospheric  effects,  and  propagation  losses  were  turned 
off.  In  case  two,  these  effects  were  turned  on.  The  radars  were  simulated  for  the  first  12  hours  after  the  collision, 
and  statistics  on  observation  errors  and  number  of  tracks  established  were  computed. 

The  radar  simulation  for  case  one  recovered  at  least  one  track  for  99  off  the  objects  in  the  first  12  hours  after  the 
collision,  with  an  average  positional  error  of  306m.  Simulation  case  two  with  all  modules  turned  on  only  recovered 
96  of  the  objects  with  a  higher  average  positional  error,  325m.  The  number  of  objects  which  had  preliminary  tracks 
determined  by  each  sensor  are  shown  in  Figure  6.  It  can  be  seen  that  the  Australian  space  fence  site  at  Canberra  was 
the  least  useful  at  quickly  recovering  tracks  based  on  this  collision. 


Objects  Tracked  at  Each  Radar  Site 
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Figure  6.  A  comparison  of  the  objects  tracked  at  each  radar  site  using  the  high  fidelity  model.  The  Gila 
River  and  Lake  Kickapoo  sites  recovered  the  most  tracks,  with  the  Lake  Kickapoo  site  recovering  the  most 
objects  not  recovered  by  any  other  sensor. 


Timing  results  were  captured  for  case  one  and  case  two.  For  case  one  using  eight  processors  per  radar,  the 
simulation  took  only  172  seconds  to  complete.  Case  two  with  all  atmospheric  effects  turned  on  took  considerably 
longer,  802  seconds  with  eight  processors  per  radar.  The  parallel  scaling  was  tested  for  case  two  from  2-16 
processors  per  radar,  with  the  timing  results  shown  in  Table  1.  It  can  be  seen  that  the  parallelization  is  around  80% 
efficient;  doubling  the  number  of  processors  causes  the  run  time  to  drop  by  40%. 


Table  1.  Simulation  Run  Time  Scaling 


Processors 

Run 

Per  Radar 

Time  (s) 

2 

2020 

4 

1265 

8 

802 

16 

487 

4.  CONCLUSION 

A  high-fidelity  parallel  radar  simulation  was  developed  for  modeling  existing  or  proposed  sensors  in  the  Space 
Surveillance  Network.  This  model  includes  antenna  effects  for  modeling  errors  due  to  the  IMU’s  and  phase  shifters, 
as  well  as  polarization  and  gain  patterns.  It  includes  antenna  temperature  noise  due  to  galactic  noise  distribution,  as 
well  as  atmospheric  effects.  These  high-fidelity  modules  are  time  consuming  to  compute,  so  the  software  is 
parallelized  on  a  per-pulse  basis.  Using  this  parallelization,  four  potential  space  fence  radars  were  tested  against  414 
debris  objects  for  12  hours  of  observation  in  less  than  10  minutes  of  wall  time.  This  speed  shows  the  potential  for 
this  radar  simulation  to  be  used  for  radar  trade  studies  and  gap  analysis  off  the  SSN  sensors. 
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